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Angular resolved electron energy-loss spectroscopy (EELS) gives access to the momentum and the energy 
dispersion of electronic excitations and allows to explore the transition from individual to collective excitations. 
Dimensionality and geometry play thereby a key role. As a prototypical example we analyze theoretically the 
case of Buckminster fullerene Ceo using ab initio calculations based on the time-dependent density-functional 
theory. Utilizing the non-negative matrix factorization method, multipole contributions to various collective 
modes are isolated, imaged in real space, and their energy and momentum dependencies are traced. A possible 
experiment is suggested to access the multipolar excitations selectively via EELS with electron vortex (twisted) 
beams. Eurthermore, we construct an accurate analytical model for the response function. Both the model and 
the ab initio cross sections are in excellent agreement with recent experimental data. 

PACS numbers: 79.20.Uv,3L15.A-,36.40.Gk 


Plasmonics, a highly active field at the intersection of 
nanophotonics, material science and nanophysics 11], has a 
long history dating back to the original work of Gustav Mie 
on light scattering from spherical colloid particles 12, 3]. For 
extended systems the plasmon response occurs at a frequency 
set by the carrier density while in a finite system topology and 
finite-size quantum effects play a key role. E.g., for a nano¬ 
shell 14-6] in addition to the volume mode, two coupled ultra¬ 
violet surface plasmons arise having significant contributions 
from higher multipoles, as demonstrated below. Such excita¬ 
tions can be accessed by optical means as well as by electron 
energy-loss spectroscopy (EELS) 17, 8]. Particle-hole (p-h) 
excitations and collective modes may “live” in overlapping 
momentum-energy domains and couple in a size-dependent 
way that cannot be understood classically 19-11]. Giant 
plasmon resonances were measured in buckminster fullerene 
C 60 112-17] and explained, e.g., by assuming C6o to have 
a constant density of electrons confined to a shell with in¬ 
ner (Ri) and outer (R 2 ) radii (the spherical shell model) 118- 
20]. Refinements in terms of a semi-classical approximation 
(SCA) incorporate the quantum-mechanical density extending 
out of the shell Ri < r < R 2 (so-called spill-out density 121]). 
Time-dependent density functional theory (TDDFT) 121-23] 
was also employed in a number of calculations 124-26], how¬ 
ever, most of them use the jellium model, i.e., the ionic struc¬ 
ture is smeared out to a uniform positive background. 

We present here, to our knowledge, the first atomistic full- 
fiedge TDDFT calculations for EELS from C 60 at finite mo¬ 
mentum transfer. We demonstrate the necessity of the full 
ab-initio approach by unraveling the nature of the various 
contributing plasmonic modes and their multipolar character. 
This is achieved by analyzing and categorizing the ah ini¬ 
tio results by means of the non-negative matrix factorization 
method 127]. The results are in line with recent experimen¬ 
tal findings 128]. The analysis also allows for constructing an 
accurate analytical model response function. 

In first Bom approximation for the triply-differential cross 
section (TDCS) for detecting an electron with momentum p/. 
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i.e., measuring its solid scattering angle dQ and energy 6p^ is 


d^cr 

dtudQ 


r Pi 


( 1 ) 


Here, p/ is the incidence momentum corresponding to an en¬ 
ergy 6pp y is the Lorenz factor, q = py^ - p/ is the momentum 
transfer, and oj = - 6p. (atomic units are used through¬ 

out). ^(q, tu) is the dynamical structure factor akin solely 
to the target 129]. The fluctuation-dissipation 129] theorem 
links S (q, of) with the non-local, retarded density-density lin¬ 
ear response function ;^^(r, r'; r - t') 129-31] via S(q,oj) = 
-(l/7r)Im[Y^(q, -q; oj)]. On the other hand,;if^(r, r'; t-f) de¬ 
scribes the change in the system density Sn(r, t) upon a small 


(a) (b) 



FIG. 1. (Color online) The T’-resolved constituents of the dynamical 
structure factor of C 60 , \lYa[6nfiq, cu)]! for (a) ^ = 0, (b) ^ = 1, (c) ^ = 
2, (d) 7’ = 3, and (e) I - A. The C 60 molecule was treated in standard 
truncated icosahedric geometry with bond lengths rc-c = 1.445 A 
and rc=c = 1.390 A. 
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perturbing potential ^^(r, t), i.e. 


dn{Y, t) ■■ 




(r, r'; t - t')S(f(r', f). 


( 2 ) 


The response function is determined by evaluating the den¬ 
sity variation with tunable perturbations, as accomplished 
via TDDFT which delivers dn{Y,t) upon solving the time- 
dependent Kohn-Sham (KS) equations [32]. 

Along this line, we utilized the Octopus package [33, 34], 
and propagated the KS equations. Kohn-Sham states are rep¬ 
resented on a uniform real space grid [35] (0.2 A grid spac¬ 
ing) confined to a sphere with 10 A radius. For the ground 
state we checked the performance of different typical func¬ 
tionals and found that the local-density approximation (LDA) 
improved by self-interaction correction (SIC) yields fairly 
good results. The HOMO (-9.2 eV) is located slightly too 
low with respect to the experimental value (-7.6 eV) [15]. 
The band width (which is typically underestimated in DFT) 
within the LDA-rSIC scheme is the largest for the tested 
functionals [36]. LDA-type Troullier-Martins pseudopoten¬ 
tials are used to incorporate the influence of the two core 
electrons per C atom, such that only the 240 valence elec¬ 
trons accounted for. Gaussian smearing has been employed 
to deal with the degeneracy of the HOMO. In gas-phase the 
molecules are randomly oriented. Hence, we have to evalu¬ 
ate the spherically averaged structure factor S{q, co). Techni¬ 
cally, this can be accomplished by choosing the perturbation 
d^(r, t) = IoS(t)j^(qr)Y*^(^r) [37] where is the spherical 
Bessel function and is the spherical harmonic. The pertur¬ 
bation strength lies with Iq = 0.01 a.u. well within the regime 
of linear response. The perturbed states are then propagated 
by the AETRS propagator [38] up to T = 20 fi/eV with a time 
step of Ar = 2 X 10“^/z/eV, covering the range from 0.31 eV to 
3142 eV in frequency space. The large simulation box ensured 
the adequate representation of excited states. A mask was 
multiplied to the Kohn-Sham states at each time step in order 
smoothly absorb contributions above the ionization threshold. 
From the density variation Sn(r,t) = n(r,t) - n(r,t = 0), 
= f drd^(r, 077(^^)f7m(f^r) is then computed in 
each time step and Fourier transformed to Sn^miq, oj) allow¬ 
ing to determine S (q, oj) as 


^ ^max ^ 

S{q,Oj) = 2j w)] 

0 f—f\ p 


( 3 ) 


£=0 m=-£ 


The m-dependence is subsidiary. To a good approximation 
henceforth m = 0 (cf. Eq. (3)). It is sufficient to consider 
\\m{dni{q, (jS)\\ = -Im[d^^<^)] which stands for the i- 
resolved dynamical structure factor depicted in Eig. 1. Eor 
^ ^ 0 (in the optical limit) the dipolar term is clearly domi¬ 
nant over higher multipoles. 

According to the shell model [20] the C6o molecule possesses 
a volume plasmon mode (^ = 0 and radial density oscillation 
with one node), a symmetric surface mode {i > \ and no ra¬ 
dial oscillation), and an anti-symmetric surface mode {i > \ 
and one radial node). We denote these modes by V, SI and 
S2, respectively. The plasmon energies are derived as coy = 

^(1,2)/ = #[1 + 2^V1+W+1)(/?i/«2F^‘>]. 



FIG. 2. (Color online) (a) Frequency-dependent part of the SI modes 
obtained from the NMF (shaded curves) with fits (dashed lines). For 
plasmon features we concentrate on the region cjj > 18 eV. (b) q- 
dependent part of the SI modes from the NMF (solid lines) along 
with fits using the model fluctuation density (symbols), (c) Model 
fluctuation density psi,^ in a plane cut through the center of the 
molecule for i ranging from 1 (top) to 4 (bottom). 


Inspecting the 7” = 1 panel the two surface modes may be 
identified around q ~ 0.3 A"\ u; ~ 20 eV and q - \ A"\ 
u; ~ 40 eV. 

As evident from Fig. 1, for higher q plasmonic modes 
(SI, S2, V) seem to merge and attain various multipoles con¬ 
tributions. This is a manifestation of electronic transitions be¬ 
tween the single-particle states with different angular momen¬ 
tum [39-41]. Thus, the question arises of how to disentangle 
these modes and to unravel their multipolar nature. 

A suitable mathematical tool to tackle this task is the 
non-negative matrix factorization (NMF), which is exten¬ 
sively used, e. g., for face recognition algorithms [27]. 
Applied to our problem, the NMF delivers two functions 
Fi{cjS) > 0 and Gi{q) > 0 that enter the density response 
as \Sn^{q, oj)\ = 2/ Fi((jj)Gi(q) (see appendix A). This struc¬ 
ture follows namely from the Lehmann representation of 
X^(,Y,r’\a)) as 


r'; w) = y ^„(w)p„(r)p<,(r'), —— 

^ (w + iTaf - El 

( 4 ) 

where pa is the real fluctuation density corresponding to a 
transition from the ground to an excited many-body state la¬ 
belled by a (with excitation energy Ea), and is the line 
width. Assuming spherical symmetry, excitations have angu¬ 
lar (7) and radial (v) components. Expanding pa(Y)pa(Y') = 
'LemRv,e{r)Rv,e{r'Wemi!^v)Y*(„{^v') Eq. (4) implies for the 
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FIG. 3. (Color online) (a) Frequency-dependent part of the S2 modes 
from the NMF (shaded curves) and corresponding fits (dashed lines). 
For the latter, no constraint has been imposed on the frequency range, 
(b) ^-dependent part of the S2 modes from the NMF (solid lines) and 
fits (symbols), (c) Model fiuctuation density ps 2 / as in Fig. 2. 


structure factor 

S(q,oj) = + l)FyA0J)GyAq), 

Gy/(q) = I^J' drr^RyAr)M‘ir)j ■ 

In full generality the sum (4) contains infinite number of terms 
corresponding to the infinite number of excited states. For ho¬ 
mogeneous electron gas plasmons are strongly damped when 
their momentum enters the p-h continuum, where the non¬ 
interacting structure factor S^^\q,CL>) > 0. For electrons con¬ 
fined to a spherical shell the momentum can be represented 
by a magnitude q and an angular momentum To mark the 
effective region and ^max in which plasmon modes exist, 
we estimate the transverse momentum as linlR (with radius 
R) and compare it to the critical momentum ^crit = 0.559kF 
[42] (the Fermi momentum is kp = ( 97 r/ 4 )^/'^r“^). We find 
so a critical ^ ~ 3. Thus, any collective excitation beyond 
^max = 4 will be suppressed. For a complementary picture, 
we analyzed S^^\q,(jS) in SC A [43], for which the electron 
density enters as a central ingredient (we take the spherically- 
averaged DFT density ^o(^)) [44]. This allows to estimate for 
which q the p-h pairs dominate the spectrum for each £ sepa¬ 
rately. For ^max = 4 we find the p-h domain at ^ > 1.2 A“^ 
Note that due to geometrical confinement plasmons and p-h 
excitations intersect each other and couple so significantly. 
Now we separate the response into v = SI (Fig. 2) and 
y = S2 (Fig. 3) for ^ > 1, while the mode y = V can 
be found from ^ = 0 density component (Fig. 4). The on¬ 
set of p-h excitations is also present in the spectra. The 
plasmon frequencies ojy^i are identified from the maximum 



FIG. 4. (Color online) (a) Frequency-dependent part of the VI (blue 
shaded curve) and V2 (purple shaded curve) mode from the NMF 
along with corresponding fits (dashed lines). Fitting has been carried 
out in the complete frequency range, (b) ^-dependent part of the V1 
and V2 modes from the NMF (solid lines) and fits (symbols), (c) 
Model fluctuation densities in the same plane as in Fig. 2. 


of the tu-dependence spectra as obtained by the NMF in the 
form I{{cd -r iTy,^)^ - Inspecting 

Fig. 2(a), we find the dipole plasmon at ojsi,i - 21.59 eV, 
this is a well established value. Increasing ^ shifts the peak 
to larger energies (in line with the shell model); the sharp 
peak around 7.5 eV, which is known to consist of a series 
of p-h excitations [ 11 ], gains spectral weight until it dom¬ 
inates for 7” = 4. Abundance of large angular momentum 
states around HOMO-LUMO gap [41] increases the number 
of channels for high-multipole electronic transitions and is 
responsible for the peak’s enhancement. The plasmon fre¬ 
quency 6 g)si ,4 = 25.64 eV on the other hand is smaller than 
^si ,3 = 26.03 eV. This demonstrates the limitations of the 
SCA. 

The radial profile of the density oscillations /?y/(r) can 
be inferred from Gy/(^) in that we assume Rl\^(r) = 

A^r exp[-(r - r^)^/ 2 (T^] and extract the parameters (A^, cr^) 

for which \\G^^ ^(q) - dr ^(r) jfiqr)^ || is minimized. 
The effective fluctuation densities are then given by psi/(r) = 
RsiA^Wmi^r), cf. figure 2 (c). 

Analogous procedure for S2 modes (Fig. 3) reveals a de¬ 
crease of the plasmon energies in qualitative agreement with 
Ref. [14]. However, the dispersion is less pronounced than 
in the shell model. To characterize the fiuctuation densi¬ 
ties, we use an Ansatz containing a node ^ 52 = A^r(l - 

r/rf^) exp[-(r - r^)^/ 2 (T^] and determine the parameters as to 
match Gs 2 ,^(^) (Fig- 3(b)). The spatial structure of the plas¬ 
mon oscillation is shown in Fig. 3(c). 

A common and physically intuitive feature of the SI and 
S 2 modes is that the spatial extend of the fiuctuation density 
is growing with £. This is a consequence of the increasing cen¬ 
trifugal force, ’’pushing” the oscillation away from the center. 

Applying the NMF with two components to |Im[^^o(^, ^)]| 
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FIG. 5. (Color online) TDCS for EELS of C6o at scattering angles 6 - 
3° (a), 6-4° (b), and 6 = 5° (c). The energy loss o) is with respect 
to initial beam energy of eq = 1050 eV. Colored curves represent 
TDDET calculations resolved in the contribution from SI, S1+S2 
and S1+S2+V. The thick curve shows experimental data [28]. (d-f): 
comparison of full TDDET and model cross sections. 


shows (Fig. 4) that in addition to the expected volume plasmon 
(labelled by VI) around tuyi = 42.69 eV (which agrees well 
with density parameter ~ 1), a second resonance peaked 
around tuv 2 = 24.17 eV appears. To clarify its origin we com¬ 
puted the response function from its non-interacting coun¬ 
terpart in the random-phase approximation and invoking the 
SCA (see appendix B). After obtaining |Im[^^o(^. tu)]| we ap¬ 
plied the NMF, as well. This procedure yields very similar 
spectra including the occurrence of V2. This feature is, how¬ 
ever, very sensitive to the details of the density distribution; it 
vanishes for a discontinuous step-like profile. Thus, it is the 
oscillations of the spill-out density taking place on the surface 
of the molecule that form V2. This is a pure quantum effect. 


also observed for the SI modes due to their dispersion with 
respect to 

Furthermore, the extracted tu-dependencies and the model 
fluctuation densities can be used to construct an approximate 
structure factor + l)F^^^(oj)Gf^(q) that 

reproduces the TDDFT results around the plasmon resonances 
in a precise way by construction. Corresponding TDCSs are 
compared in of Figs. 5 (d-f). 

An important feature of the structure factor is the /-sum 
rule J^dojcjijSiq^oj) = Nq^l2 (number of electrons N). 
Checking for the (plasmon-dominated) oj) shows the 

discrepancy for larger q\ a critical value is reached when 

dciJCijS^^^^\q,co) decreases again after quadratic growth. 

We find ^crit ~ 1.2 A“^ which is consistent with the estima¬ 
tion above. Hence, p-h excitations become more important 
for q > ^crit and gradually diminish the plasmon contribution. 

In summary, we presented accurate TDDFT calculations 
for the dynamical structure factor and EELS spectra for C6o 
molecule underlining the role of higher multipole contribu¬ 
tions. Using NME decomposition allowed to trace the evo¬ 
lution in q and tu of the symmetric and anti-symmetric sur¬ 
face and volume plasmons. In addition, we characterized 
and modeled the fluctuation densities (i.e., the ingredients of 
the response function) and unveiled their multipolar character. 
These ingredients might, in principle, be accessed selectively 
by using electron beams carrying a definite angular momen¬ 
tum (electron vortex beams [45, 46]). By measuring the angu¬ 
lar momentum of the scattered beam the angular momentum 
transfer M becomes a control variable which the EELS spec¬ 
tra depends on. Particularly, provided the beam axis coincides 
with the symmetry axis of spherical system, the plasmonic 
response upon scattering of such twisted electrons contains 
multipole contributions for 7” > |A^| only [47]. Hence, specific 
multipoles can be excluded or included by varying M. 

Eurthermore, we discussed the limitation of spherical-shell 
models in describing the quenching of the volume plasmon 
and identified the electronic density distribution as a key factor 
determining its energy. We obtained excellent agreement with 
experimental results and explained how the different plasmon 
modes contribute to the spectra. 


With the dynamical structure factor being fully character¬ 
ized, we proceed by computing the TDCS (Eq. (1)). Eig. 5 
compares calculated and measured [28] EELS spectra as a 
function of the electron scattering angle 6 which fixes the 
momentum transfer. The magnitudes of the measured spec¬ 
tra shown in Eig. 5 are determined up to an overfall factor. 
Thus, the theory-experiment comparison in Eigs. 5 (b,c) is on 
an absolute scale. The classification of the plasmon modes 
accomplished by the NME analysis allows for plotting mode- 
resolved TDCS curves. As Eigs. 5(a-c) demonstrate, the SI 
plasmons play the dominant role for small 6 (which corre¬ 
sponds to the optical limit of small q), while the S2 modes 
becomes increasingly significant for larger 6 (i.e., larger q). 
The larger energy of the S2 with respect to the S1 plasmons 
leads to the formation of a shoulder (clearly visible for 0 = 4°) 
and, thus, to the apparent shift of the maximum of the exper¬ 
imental EELS spectrum with growing 0. A similar effect is 


Appendix A: Non-negative matrix factorization 

As dictated by the fluctuation-dissipation theorem, the 
imaginary part of Sn^iq, oj) for tu > 0 is purely negative. Thus, 
the non-negative matrix factorization (NME) can be applied to 
\lm[Snf(q, oj)]\ = -lm[Snf(q, oj)] to split 

N 

|lm [6ne(q, w)]| = ^ Fy^e{(^)Gy/{q). (Al) 

y=l 

Without imposing any restriction on the number of compo¬ 
nents (A) the expansion (Al) is exact and can be paralleled 
with the singular value decomposition (SVD) of a general 
(complex or real) matrix M: M = UEV*. The difference is in 
the additional requirements of positivity on the vectors form¬ 
ing U and V. The transition from continuous variables as in 
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eq. (Al) to the matrix form is provided by discretizing the oj- 
and the ^-points after smooth interpolation. 

We select N = 2 as we expect two dominant surface 
plasmon modes (SI and S2). This choice is confirmed by 
computing the residue norm with respect to the full function 
\lm[Snf(q, oj)]\. 

The problem of non-negative matrix factorization can be 
formulated as a non-convex minimization problem for the 
residue norm r = IIA - WH|p. Thus, the solution is not 
unique and may lead to local minima. Depending on the norm 
used different algorithms can be formulated. A commonly 
used method is the multiplicative update of D. Lee and S. Se- 
ung [27]: 


W < W 

“ “(WHH^),Z 

(A2a) 

(W^A)„,- 

(A2b) 


where i indexes the energy points and j numbers the time 
points. The method starts with some suitable guess for ma¬ 
trices W and H. Additionally, the vectors forming W are nor¬ 
malized each step: 


u 


iiwjr 


Upon these prescriptions (A2) the Euclidean distance r 
monotonously decreases until the stationary point (local min¬ 
imum) has been reached. We initialized the vector Wi (W 2 ) 
with cuts of ^)]| along q direction at u; = 20 eV 

(oj = 40 eV), while Hi (H 2 ) is constructed by cuts at ^ = 
0.5 (q = 1.0 A“A- We found that typically 1000 itera¬ 
tions yield well converged results. 

The functions Fy^^ioj) and Gy^^iq) is then obtained from in¬ 
terpolating the data from Hy and Wy, respectively. We nor¬ 
malize the frequency spectra such that fitting by = 

Im[26t;y//((6i; -r iTy^^)^ - col^)] (as explained in the main text) 
can be performed without any additional prefactor. Gy^^iq) is 
normalized accordingly. This normalization procedure is con¬ 
sistent with the Lehmann representation. 


point is the Dyson equation for the density-density response 
function in random-phase approximation (RPA): 


Xir:,r:'-,z) = V°\r,r';z)+ 



dr2V°^(r,ri;z) 


xu(ri -r2Mr2,r';z). 


(Bl) 


We drop the superscript R and consider general complex argu¬ 
ment z here. In SC A, the non-interacting reference response 
function r'; z) can be expressed in terms of ground-state 
density no(r) (= no(r) as we assume spherical symmetry here) 
only. The subsequent derivations and the solution scheme for 
eq. (Bl) are detailed in the Supplementary Material [48]. The 
amount of spill-out density can be adjusted by varying the 
smearing parameter Ar in the model density 


noW = No[6Ar(r - Ri) - 9Ar{r - 7?2)] , 

1 (B2) 

e^rir) = —-—— , 

1 + exp[-?-/A?-] 

where Ri = Ro-AR/2, R 2 = Rq + AR/2 are the inner and outer 
radii (Rq = 6.5 a.u.), while the normalization No ensures the 
correct total valence charge. AR is fixed to keep the mean den¬ 
sity constant. The scenario Ar ^ 0 corresponds to a box-like 
density profile with sharp boundaries, while Ar = 0.5 a.u. is a 
good approximation to the spherically-averaged DFT density. 
Once eq. (B 1) is solved for certain Ar, the = 0) contribution 
to the structure factor, lm[Sno(q,z = oj iT)] (T = 0.1 a.u. is 
a broading parameter) can be computed. Applying the NMF 
technique allows again for separating the VI and V2 modes. 
We find the position of VI similar to the TDDFT results, while 
the behavior of V2 is very sensititve to Ar. While very pro¬ 
nounced for Ar = 0.5 a.u., the relative strength of the V2 peak 
vanishes for Ar ^ 0. More details and graphs of volume plas¬ 
mon spectra can be found in the Supplementary Material. 
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